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ABSTRACT 



We present results of the investigations of the statistical properties of a joint 
density and velocity divergence probability distribution function (PDF) in the mildly 
] non-linear regime. For that purpose we use both perturbation theory results, extended 

^ ■ here for a top-hat filter, and numerical simulations. 

[~ — ' In particular we derive the quantitative (complete as possible up to third order 

l/^ , terms) and qualitative predictions for constrained averages and constrained disper- 

' sions - which describe the nonlinearities and the stochasticity properties beyond the 

linear regime - and compare them against numerical simulations. We find overall a 
, good agreement for constrained averages; however, the agreement for constrained dis- 

persions is only qualitative. Scaling relations for the fi-dependence of these quantities 
are satisfactory reproduced. 

r*^ . Guided by our analytical and numerical results, we finally construct a robust 

CIh' phenomenological description of the joint PDF in a closed analytic form. The good 

agreement of our formula with results of N-body simulations for a number of cosmo- 
^ • logical parameters provides a sound validation of the presented approach. 

^ ' Our results provide a basis for a potentially powerful tool with which it is possible 

^ , to analyze galaxy survey data in order to test the gravitational instability paradigm 

beyond the linear regime and put useful constraints on cosmological parameters. In 
particular we show how the nonlinearity in the density-velocity relation can be used 
. to break the so-called fJ-bias degeneracy in cosmic density-velocity comparisons. 

^ ' 

■ Key virords: cosmology: theory - galaxies: clustering - galaxies: formation - large- 

scale structure of the Universe 



1 INTRODUCTION 

Comparisons between cosmic density and peculiar velocity 
fields of galaxies have become already a small industry in 
cosmology (e.g., Strauss & Davis 1988; Yahil 1988; Kaiser et 
al. 1991; Dekel et al. 1993; Hudson 1993; Hudson et al. 1994; 
Davis, Nusser & Willick 1996; Riess et al. 1997; Willick et 
al. 1997; da Costa et al. 1998; Sigad et al. 1998; Willick & 
Strauss 1998). The goal of this pursuit is to test the gravi- 
tational instability hypothesis for the formation of structure 
in the Universe and to measure the cosmological parameter 
Q. In the paradigm of gravitational instability the density 
and the velocity fields are tightly related and the relation 
between them depends on Q. In the linear regime, i.e., when 
the density fluctuations are significantly smaller than unity, 
this relation reads. 



5{x)^-fin,Ar'\/-vix). (1) 

Here, 5 is the mass density fluctuation fleld, v is the veloc- 
ity field, f{0, A) ~ fi"'^ and we express distances in units 
of km s~^. However, the derived amplitude of the density 
fiuctuations from current redshift surveys (e.g.. Fisher et 
al. 1995) and from the potent reconstruction of density 
fields (Dekel et al. 1990 and Bertschinger et al. 1990), goes 
somewhat beyond the linear regime. For example, the den- 
sity contrast in regions like the Great Attractor or Perseus- 
Pisces is about unity even when smoothed over scales of 
1200 km s"^ (Sigad et al. 1998). Future redshift surveys 
and peculiar velocity catalogs are expected to provide reli- 
able estimates of density and velocity fields on scales where 
nonlinear effects are certainly non-negligible and may lead 
to interesting consequences such as breaking the degener- 
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acy between and bias (Chodorowski & Lokas 1997, here- 
after CL97; this paper; Chodorowski, in preparation) and 
therefore need to be accounted for in performed analyses. 
To date, there have been several attempts to construct a 
mildly nonlinear extension of relation (jl|) . Those were either 
based on various analytical approximations to nonlinear dy- 
namics (Nusser et al. 1991; Bernardeau 1992b; Gramann 
1993; MancineUi & Yahil 1995; Chodorowski 1997), or N- 
body simulations (MancineUi et al. 1994; Canon et al., in 
preparation) . 

The aim of this paper is to describe quantitatively the 
the density versus velocity-divergence relation (DVDR) at 
mildly nonlinear scales.^ We focus on that regime because 
it is explored in current analyses of the observational data, 
which are commonly smoothed over scales chosen to ensure 
that the r.m.s. density fluctuation, though significant, is be- 
low unity. Also this is the domain of applicability of PT, 
permitting us applying some of the recent results obtained 
within its framework and relevant to our goal. 

Recently, CL97 and Chodorowski et al. (1998, hereafter 
CLPN) derived the DVDR up to third order in perturbation 
theory (PT), assuming Gaussian initial conditions. Specif- 
ically, CL97 derived the 'forward' relation (the density in 
terms of the velocity divergence), while CLPN derived the 
'inverse' relation (the velocity divergence in terms of the 
density) and the scatter in both relations. Due to the scat- 
ter, the inverse relation is not identical with a mathematical 
inversion of the forward one. The coefficients in these rela- 
tions were calculated for the case of Gaussian smoothing of 
the evolved fields which is also applied in the analysis of ob- 
servational data. For large smoothing scales, N-body results 
are expected to converge to the results of CL97 and CLPN. 
Indeed, the predicted values of the coefficients proved to be 
in qualitative agreement with the results of Canon et al. (in 
preparation), who tested similar functional forms of the re- 
lations against N-body simulations. On the other hand, for 
small smoothing scales (at which the r.m.s. density fluctua- 
tion is close to unity), contributions from orders higher than 
third may become significant. 

In this paper, instead of developing a PT methodology 
beyond the third order, we make use of N-body simulations. 
Yet, in the mildly nonlinear regime, our approach is still 
superior to those based solely on means of N-body simula- 
tions. The asymptotically exact formula, we wish to modify 
here is the complete third-order relation of CL97 and CLPN, 
not the linear expression given by Equation (|l|). Indeed, as 
we will see the modification to the third-order results con- 
cerning the mean (the 'forward' and 'inverse') relations is 
quite modest. In other words, the third order mean relations 
prove to be already good approximations. However, since 
the mildly nonlinear DVDR is not deterministic, its full de- 
scription should include not only the constrained averages 
but the whole joint Probability Distribution Function (PDF) 
for density and velocity divergence. We construct here this 
PDF in a closed analytic form. 



* We define a density field to be mildly nonUnear if the r.m.s. 
density fluctuation is a significant fraction of, but still smaller 
than, unity. Then the mildly nonlinear scales in the Universe are 
these about or greater than 8 Mpc. 



In order to establish accurately the DVDR from N- 
body simulations, one has to treat properly the final veloc- 
ity field, determined in a simulation only at a set of discrete 
points (final positions). Commonly used smoothing proce- 
dures lead to effectively a mass averaged velocity field, which 
induces spurious velocity gradients. The problem of proper 
volume-weighting of the velocity field was solved for a top- 
hat smoothing by Bernardeau & van de Weygaert (1996) 
using Voronoi and Delaunay tessellations. Therefore, here 
we will investigate the DVDR for top-hat smoothed fields. 
The results of CL97 and CLPN were obtained for Gaus- 
sian smoothing; consequently we extend them for top-hat 
smoothing. The algorithm of Bernardeau & van de Wey- 
gaert (1996) can be possibly modified so as to apply it to 
Gaussian-smoothed fields. This will be addressed elsewhere 
(Chodorowski & Stompor, in preparation). 

In sum, new parts of this investigation are: 

(i) Calculation up to third order in PT of the coefficients 
of the mean relations and of the scatter in both, for top-hat 
smoothing; 

(ii) Comparison of PT predictions to N-body simulations 
in the mildly nonlinear regime; 

(iii) Construction of the joint PDF for density and veloc- 
ity divergence. 

Also, we will show how the nonlinearity in the DVDR can be 
used to break the so-called fi-bias degeneracy in comparisons 
between large-scale density fields of galaxies and the fields 
of their peculiar velocities. 

In Sect. ^ we present the results obtained from PT for 
the constrained averages. Sect. ^ is devoted to the presenta- 
tion of numerical results. We first confront them against PT- 
derived predictions and subsequently guided by both kind 
of results construct a phenomenological description of the 
joint density-velocity PDF. The technical calculations are 
deferred to appendices. These contain in particular the ex- 
plicit expressions of the coefficients that intervene in PT. 



2 PERTURBATIVE RESULTS 
2.1 General features 

At large scale, it is expected from linear theory that the 
local density contrast and the velocity divergence are pro- 
portional (Eq. 0). However, as soon as nonlinear effects come 
into play, and this is expected to happen when the variance 
reaches a fraction of unity, this relation is expected to fail. 
The true relation departs from that given by Equation (^ 
in a immber of ways. First, its linearity evidently does not 
hold any more (see Fig. |l]) . Secondly, and more importantly, 
no more a one-to-one correspondence exists between local 
density and local velocity divergence. Non-local effects in- 
troduce a scatter in this relation. As a result, the descrip- 
tion of the relation between the local density and the local 
convergence should be made with the construction of the 
joint Probability Distribution Function (PDF) of these two 
quantities. Let us define the velocity divergence field, 

e{x) -vix), (2) 

and the scaled velocity divergence with 

e(x) = ~f(n,A)-' e. (3) 
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Figure 1. A typical joint PDF obtained in numerical simulations 
(see text). The smoothing scale is 15/i~^Mpc, and the contour 
plot is made in logarithmic scale. 



Note the minus sign and the factor /(n^A)"^ muhiplying 
the velocity divergence. Linear relation nw is then simply 
5 — 9. In the linear regime this PDF, Pjoint, therefore reads, 



P/Zr(5,e) dS de^exp 



2a2 



r ,r ~, d5 d6' 

foirac(<?-6>) (4) 

V 2 TT (T^ 



where a is the r.m.s. fluctuation of the density field. Non- 
linear couplings will not only introduce some non-Gaussian 
features in the PDF of each of the variables, they will also 
change the Dirac delta function ((Joirac) in a more compli- 
cated function of a finite width. As long as the departure 
from the linear regime is small, generic features can be in- 
ferred from Perturbation Theory (Bernardeau 1992b). The 
r.m.s. value of the scatter around the mean relation between 
5 and 6 is expected to be of the order of a^. (As stressed by 
Chodorowski 1997, this can be interpreted as an effect due 
to the couplings of the local density with the shear field.) 
Therefore, for a beyond linear regime but still below unity 
the most probable values of {6, 9) form an elongated region 
in a (5, 9) plane as it is evident in Fig. ^. 

The mean trend of the DVDR is described by con- 
strained averages, either by the mean density given 9 or, 
equivalently, the mean 9 given 5. The dispersion around that 
mean trend provides us with the complementary information 
that is the amplitude of the scatter. All these features are 
a priori accessible to analytic calculations in PT; some of 
them have been done in the previous papers (Bernardeau 
1992b, CL97, CLPN). The results however are dependent 
on the chosen window function. Since filtering in Voronoi 
tessellations has only been implemented for top-hat window 
function we decide to focus our presentation on this case. 



2.2 Mean relations 

The generic formal expression of the mean relation is a stan- 
dard calculation in statistics. We sketch its derivation in 



Appendix^. In general the constrained mean density, 

is a function of both 9 and a. These two variables are a 
priori small and of the same order, however to allow us the 
description of rare events it is convenient to separate the ex- 
pansion in 9 from the one in a. An important result obtained 
by Bernardeau (1992b) is the expression of the constrained 
density as a function of S in a case of a vanishing variance. 
The constrained means have been found to be given by the 
5-9 relation exhibited in the spherical collapse model. A sim- 
plified useful expression (which is strictly valid in the limit 
f2 — > 0) is given by. 



1 + y 



3/2 



- 1, 



for 



0. 



(5) 



Though this result was obtained without taking into 
account the filtering effect, progress made in understanding 
of filtering has shown that in a case of a top-hat filter the 
filtering effects can be described as a simple mapping from 
the Lagrangian to the Eulerian space (Bernardeau 1994b). 
As a consequence the relation is expeeted to he valid for 
hop-hat filtering as well. 

When the limit cr —> is dropped only a finite number 
of terms in a joint expansion in a and 9 is known. They have 
been obtained in several recent papers (Bernardeau 1992b, 
CL97, CLPN). Up to the third order the mean S given 9 can 
be written generically as. 



((5)|- = ai{ag) 9 + a2(P - cr|) + 036 
and mean 9 given 5, 



{0)\ 



ri(as)5 + r2{S ~ as) + r^S 



(6) 



(7) 



Here, ct| and (j| are the variances of the density and the 
scaled velocity divergence field respectively and these 
terms are included so that the global averages of S and 9 
vanish. To have a consistent expansion up to the third order 
both oi and ri have to be calculated up to order whereas 
the other coefficients can be computed at their leading order. 
The coefficient 02 was computed by Bernardeau (1992b) for 
a top-hat filter and it is easy to check that r2 — —02. More 
generally the coefficients Um and were derived by CL97 
and CLPN respectively, and their values were explicitly cal- 
culated for a Gaussian window function. Up to the third 
order, the coefficients a2, 03, r2 and ra do not depend on 
the normalization of the power spectrum but they depend 
slightly on the smoothing scale through the power spectrum 
index (unless it is a power law). In Appendix^ we explicitly 
compute the values of the coefficients for a top-hat window 
function. Unlike for a Gaussian smoothing, their values turn 
out additionally not to depend on the spectral index, as ex- 
pected if the picture of the Lagrangian to Eulerian mapping 
is correct. Specifically we have. 



a2 = — -0.190, 

and 



r2 



4 
21 



-0.190, 



as = 



40 
3969 



328 
3969 



-0.0101 



0.0826 . 



(8) 



(9) 



It is worth noting that 02 and are very close to those 
obtained Taylor expanding eqn. (p|) in terms of 9, indicat- 
ing extremely weak fl-dependence of the relation between 
density and scaled velocity. The values of the coefficients ai 
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and ri are at the leading order both equal to unity. Their 
next-to-leading-order corrections, as predicted by PT, are 
propo rtion al to t he va riance of the respective fields [equa- 
tions (Bff) and ( [Bf7| )]. fn fact, up to the third order it is 
sufficient to consider the linear variance of the fields, which 
is identical, and will be denoted a^. The coefflcients of pro- 
portionality of the corrective terms depend on the spectral 
index, similarly to the case of a Gaussian smoothing. We list 
the values of (ai — l)/cr^ and (ri — l)/cr^ as functions of the 



spectral index in Table CI 



2.3 Scatter in the relations 

The mean relations, both forward and inverse, have a scat- 
ter. The r.m.s. value of the scatter is equal to the square 
root of the conditional variance, 

^sio^{iS-{3)\s)X' (forward), 
^0i6 = {iS-{S)\s)X' (inverse). 

Hence up to the next-to-leading order (CLPN), 
'^419 = ^0<t| + bial§ + b2alP + 0{al) (10) 
and 

= Soaj + SiajS + S2Cr|(5^ + ©(erf) . (11) 
PT predicts furthermore that 

So = bo and S2 = ^2 ■ (12) 



Explicit formulas for the coefficients bo and &2 were derived 
by CLPN. In expression ([ic|), the terms boO"| and b2(Jg9^ 
are of the leading, 0{ai), order (as 6 ~ ag), while the term 
bia^O is of the next-to-leading order, 0((t|). We have in- 
cluded this term because for a top-hat filter the value of the 
coefficient 62 is exactly zeroQ (see Appendix^). Therefore, 
the lowest-order non-zero contribution to the ^-dependence 
of the constrained variance is already provided by a next-to- 
leading-order term. In Appendix ^ we calculate the values 
of the coefficient 60 for power-law spectra. The results are 
presented in the last column of Table CI. We found that 
the expression for 61 (and si) contains perturbative contri- 
butions up to the third-order. However, the mixed moments 
which enter the calculations (with terms of the form {V1V2V3) 
where v stands for either S or 6) are of extra complexity in 
comparison to any other moments computed here. In the 
present paper we will not therefore try to predict the value 
of 61 or si. 



3 NUMERICAL RESULTS 

The aim of this section is twofold. We want 

• to check, qualitatively and quantitatively, the PT re- 
sults presented in the previous Section. 

• to build a complete phenomenological description of the 
joint PDF of the local density and velocity divergence that 
would be valid for typical events, S of order ag and 6 of order 
ag, as well as for the rare event tails. 

t In the case of a Gaussian smoothing, 62 is not exactly zero but 
still very small. 



Except for the formula (|5|), all perturbative results were 
derived for 5 ~ erf and 6 ^ ag. For applications of these re- 
sults to current density-velocity comparisons this assump- 
tion is justified. The volume of current peculiar velocity sur- 
veys is too small, sampling too sparse and noise too big to 
detect unambiguously 'rare events' in the observed velocity 
field. However, in current N-body numerical simulations the 
simulation volume can be made significantly greater, sam- 
pling is dense and there is essentially no noise in the velocity 
field traced by particles. Therefore, the use N-body simula- 
tions enables one to leave the regime of 'typical events' in 
the (5, 6) plane, i.e., to study 'rare events' as well. This is not 
only of academic value, since with accumulation of new ve- 
locity data and discovery of new distance indicators we may 
hope the observed velocity field to improve qualitatively in 
the future. 

Here we use a set of numerical N-body simulations 
to study the DVDR in the mildly nonlinear regime (i.e., 
as < 1), but including also its behaviour in the low- and 
high-density tail. This will be especially important for con- 
structing the joint PDF for the density and the velocity 
divergence. A brief description of the N-body simulations as 
well as numerical al gori thms employed in their analysis are 



given in Subsection 3.1 



In Subsection 3.2 



and 3.3 we study 



the constrained averages and the constrained dispersions, re- 
spectively, and compare them to the PT predictions. Subsec- 
tion 3.4 provides details of how we build a phenomenological 



description of the joint PDF using both PT and simulation 
derived information. 



3.1 The simulations and the numerical algorithms 

The N-body simulations were carried out using a grav- 
ity solver based on AP^M code developed by Couchman 
(Couchman 1991), which we adapted for use in cases with 
arbitrary negative curvature (open models) and/or cosmo- 
logical constant (see for example Peebles 1980 for general 
formulas and Peacock & Dodds 1996 for a description of 
necessary changes). 

In each case, we used 128'^ particles within the cubic 
200/i"^Mpc size box with a 128x 128x 128 FFT (Fast Fourier 
Transform) grid. Both amplitudes and phases of an initial 
density field were randomly drawn from the Rayleigh and 
flat distribution respectively and initial positions and ve- 
locities of the particles were fixed using standard Zel'dovich 
movers accordingly. The starting time of the simulations was 
chosen to ensure that an average initial particle displacement 
was less than one third of the elemental cell size with the 
power on that scale being initially of the order of 2-^5 x 10~^ . 
The gravitational force smoothing was always set equal to 
the single cell size at the beginning of the simulation, kept 
subsequently constant in physical units to become eventu- 
ally frozen on one tenth of the comoving cell size ever since 
reaching that limit (Couchman 1991). 

The equations of motion were then integrated in 500 
time-steps over an expansion factor of 30 and ended at the 
time when a (ensemble averaged) standard deviation of lin- 
ear mass perturbations reached a required value at a given 
scale (usually unity on the scale of 8/i~^Mpc). 

We had at our disposal 5 simulations done for power 
law spectra of spectral index n = —1.5 and corresponding 
to different cosmological parameters, either for an Einstein- 
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de Sitter Universe, or with D, = 0.3 and A = 0, or with 
SI — 0.3 and A = 0.7 at the final stage of the simulation. We 
choose index n = —1.5 so that the results of PT that involve 



loop terms are finite (e.g. Table SI). See discussions in the 
following. 

In each simulation we applied a top-hat filter of different 
radius. In each case the filtered quantities were obtained on 
50'^ grid points. The local filtered velocity divergence was 
obtained using the Voronoi tessellation technique developed 
by Bernardeau & van de Weygaert (1996). We have checked 
that the number of tracers we have used (about 60,000) was 
sufficient to provide us with robust results for the scales of 
interest. 



3.2 The constrained averages 

Fig. |l| shows a typical result for the joint PDF. From such 
data it is easy to compute the constrained averages. We then 
use the formulas (^) and (0) as parametric functions leaving 
the three coefficients and the variance as free parameters to 
fit the measured expectation values. The fitting procedure 
uses a weighting function centered on either 5 = or 6^ = 
with a width of the order of a. We have checked that our 
fitting method provided us with robust results for the co- 
efficients. To illustrate the effect of cosmic variance on the 
values of the estimated parameters, instead of plotting their 
values averaged over 5 simulations, we decided to make com- 
bined plots of the coefficients derived from each simulation. 

In Fig. ^ we present the resulting values of the linear 
coefficient between the local density contrast and the mi- 
nus divergence (in this case, not the scaled divergence). It 
explicitly exhibits the f{Q,,A) dependence that we expect 
for different cosmologies. Fig. ^ shows the fitted values of 
12, r2, as and rs for different cosmologies and smoothing 
scales. In this case we present the results for the fitting co- 
efficients between 5 and the scaled divergence 6. Note that 
with these variables we expect the coefficients a„ and r„ to 
be independent of a cosmological model. This property is in- 
deed observed in the results. Solid horizontal lines in Figs. ^ 
and ^ represent the PT prediction for the coefficients cal- 
culated at leading order. We can see that for small values 
of the variance the coefficients derived from the simulations 
converge to their leading-order PT values. For large values 
of the variance, the coefficients slightly, but systematically 
deviate from these values. 

The next-to-leading-order corrections to the coefficients 
ifc (or Vk) are given by so-called one- loop terms involving 
perturbative contributions of {k + 2)-th order. Therefore, 
here we were able to compute a one-loop correction to the 
linear value of the coefficients ai and ri only (since we per- 
formed a complete third-order calculation.). The results are 
plotted as thin dot-dashed lines in Fig. ^ and are rather puz- 
zling. While PT predicts a significant departure of ai and 
ri from their linear values, the coefficients derived from sim- 
ulations remain in fact essentially constant as a function of 
a. It does not seem to be a numerical problem. This fail- 
ure of next-to-leading order predictions is possibly due to 
the existence of divergences in the PT calculation. This is a 
difficulty that has been encountered in other one-loop cal- 
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Figure 2. The best-fit values of the coefficients ai (triangles) 
and ri (diamonds) as determined for a set of different simulations 
(interpolated by thick solid lines for Einstein-de Sitter, dotted 
lines for open Universes with A = 0, dashed lines for non-zero 
A Universes). In this case and for this plot only the displayed 
coefficients are those describing the relation between 5 and the 
(minus) true divergence, similarly to Eqs. (jij to make their 
Q and A dependence apparent. Note a very weak dependence on 
A, as expected. The thin horizontal lines are the linear theory 
predictions. The thin dot-dashed lines are the PT predictions at 
the one loop order (see Table CI). 



culations in PT (see Jain & Bertschinger 1994, Scoccimarro 
& Frieman 1996, and this paper) 

The cosmic variance is generally quite small. Not sur- 
prisingly, the effect of the cosmic variance is the strongest 
for the coefficients as and 7-3, and on the largest scales (the 
smallest a). 



3.3 The constrained dispersions 

From third order PT we expect that the constrained vari- 
ance of S given 8 can be written as a first order polynomial 
in 6 times cr|, see equation ( |lo| ) with 62 = (for details see 
Appendix ^). This equation assumes however that ag, 
the assumption which we have relaxed in this Section. When 
6 is arbitrary, but ag is smaller than unity, we can still use 
all-order PT to predict the qualitative behaviour of the con- 
strained variance. Then it turns out that at leading order 

t However, for some of us it still remains a puzzle, since we have 
at least a partial physical explanation for the presence of the term 
oc (T?S in the forward relation. In the Zel'dovich approximation, 
the Eulerian density is a third-order form in (first-order) velocity 
derivatives, and one of the third order terms can be cast to the 
form oc S^6, where is the shear scalar (Chodorowski 1998). 
One of the contributions to is then proportional to the con- 

strained average (S^9)|^, which, since at leading order the shear 
is independent from the divergence, is (2/3)cr?9 ! 
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Figure 3. The best-fit values of the coefiicients a2 and r2 (top 
panel) and 03 and ra (bottom panel) as determined for different 
simulations (conventions are as in Fig. 2). The presented here 
coefficients are those at the quadratic and cubic terms in the S—6 
relation. 

in (jg, the formula for the 'forward' variance is given by cr| 
times a power series in 9. An object of more practical in- 
terest is the square root of the variance, or the dispersion. 
We have found that a low-order polynomial provides a sat- 
isfactory fit to the dispersion and finally we have chosen a 
second-order polynomial. Specifically, we have found that 

9 



; 0.45 



(1 + .-+^.-) 



(13) 



is an accurate fitting formula for all cosmological models in 
our range of parameters. As for the dispersion of the 'inverse' 
relation, we have found that 



(^e\s ~ 0.45 (1 + 0.55 - 0.05 5^) a| 



(14) 

was a reliable fitting form. For reasons that are not clearly 
elucidated, this form is more robust (against variation of the 
cosmological parameters), when written as a function of ag 
instead of as- 

In Fig. ^ the fitting formulas (dashed lines) are tested 
against the numerical results (solid lines) and are seen to 
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Figure 4. The contour plots show the scatter of the density as 
a function of the (minus) local divergence (top panel) and of the 
local divergence as a function of the density (bottom panel). Thick 
(jagged) solid lines show the dg^g (top panel) and crg^g (bottom 
panel) obtained from numerical simulations and compared against 
analytical fits (dashed lines) proposed in the text (see Eqn. ( [l^ 
and (p^) Note the nearly symmetric scatter around the mean 
in the upper panel which is actually roughly compatible with a 
Gaussian distribution. The plots correspond to an Einstein-de 
Sitter Universe and a smoothing scales of 15fe~^Mpc. 



be very accurate. Note that the expression (^|) explicitly 
predicts the scatter in the forward relation to vanish for 
9 —> —3/2, that is when 9 gets close to its lower bound, 
although this was not required a priori in the fitting proce- 
dure. This is an interesting result since it naturally makes 
the width of the constrained distribution of 5 small when its 
constrained average gets close to —1. (Though we would ex- 
pect to obtain a similar property for the inverse relation, the 
second-order expression (g^ fails to reproduce it. This ex- 
pression however will not appear in an analytic formula for 
the joint PDF.) Note also that the fits recover the scaling of 
the dispersions with ag predicted by PT. Finally, the values 
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of the coefficients of the constant terms in expressions (jl3|) 
and are identical, as predicted by PT. On the other 
hand, they are not equal to the PT value. Comparing equa- 
tion with equation ( p^ we see that the value of the 
parameter bo derived from N-body is given by (0.45)^ ~ 0.2, 
while, for the spectral index n = —1.5 that was used, PT 



predicts about 0.05 (see Table CI). Once again, the PT re- 



sults that involve loop terms fail. 

3.4 A description of the joint PDF 

Our goal in this paper is to construct a generalization of 
Equation (^) applicable in the mildly non-linear regime. A 
natural approach would be to rebuild the joint PDF from the 
knowledge of the various moments. The mathematical con- 
text that allows to do such a reconstruction is known. It is 
called the Edgeworth expansion (see Bernardeau & Kofman 
1995 for a derivation of the Edgeworth expansion from the 
cumulants and Appendix Q . As is stressed in Appendix Q 
in our case this expansion (or any of a similar kind) is un- 
fortunately inapplicable. The reason is that at leading order 
the cross-correlation matrix of 5 and 9 is singular. Therefore 
there exists a linear combination of both variables (namely 
S — 9) which has a vanishing variance and the one-point PDF 
of which is not necessarily close to a Gaussian distribution 
even for arbitrarily small a (although it converges to a Dirac 
delta function when a —> 0). 

The construction of the joint PDF must then rely on 
further assumptions. We assume here that the conditional 
distribution function of the density for any (fixed) value of 6 
is Gaussian. This assumption is supported by the properties 
observed in Fig.^ It appears that a similar assumption for 
the reverse relation would have been less accurate. The con- 
struction of the whole joint PDF can then be done from the 
knowledge of the constrained mean, the constrained scatter 
and the one-point PDF of 9. For the later problem we took 
advantage of the existence of a very accurate analytical fit^ 
proposed by Bernardeau (1994b), 



Pom9 = 



(27r)i/2 ag 



exp 



with 



2Aa? 



d9 



A = l + y, 



(15) 



(16) 



This distribution has been found to be very accurate in de- 
scribing the PDF for the Einstein-de Sitter Universe (see 
Bernardeau & van de Weygaert 1996) and also for other 
Friedman-Robertson- Walker backgrounds (Bernardeau et 
al. 1997) if the scaled divergence is used. 

The proposed expression for the joint PDF is then given 

by 



P{S, 9) d9 dS = exp 



d9d5 



(17) 



8 This PDF has been derived strictly speaking for a n = — 1 
power spectrum. 
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-1.0 -0.5 0.0 0.5 1.0 1.5 

local scaled divergence 



2.0 



S. 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 
local scaled divergence 



S. : 




-1 \..'.' 



-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 
local scaled divergence 

Figure 5. Comparison of the isodensity contours of the joint PDF 
obtained numerically (thin solid lines) and using the phenomeno- 
logical formula (thick solid lines) proposed in the text (Eqn. (p^)). 
Thicker lines are the constrained mean of 5 from numerical results 
(solid jagged line) and from our fit (dashed line). The three panels 
correspond to three different smoothing scales, lOh^'^Mpc (top), 
15/i~^Mpc (middle) and 25 h^^Mpc (bottom) for the Einstein-de 
Sitter Universe. 
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2 - 



a - 




-1.0 -0.5 0.0 0.5 1.0 1.5 

local scaled divergence 



2.0 



S. 




-1.0 -0.5 0.0 0.5 1.0 1.5 
local scaled divergence 



2.0 



Figure 6. Same as figure g but for other cosmologies: an open 
universe with £7 = 0.3 and A = (top panel), and a zero curvature 
Universe with Q = 0.3 and A = 0.7 (bottom panel). A smoothing 
scale is equal to 15/i~^Mpc. 



As a formula for the constrained variance a^^g we adopt 
the expression As a formula for the constrained mean 

we might adopt a third-order pol ynomial in with the co- 
efflcients fitted in Subsection 3.2. However, the coefficients 



depend weakly on the smoothing scale (equivalently, on ag) 
so we should then model their (Tg-dependence. Furthermore, 
and more importantly, as we stated earlier, in this Section 
we study the DVDR including its behaviour in the tails of 
the joint PDF. A third-order polynomial was overall a good 
fit to the data, but it did not satisfy a specific constraint 
in the negative tail. Namely, for — > —3/2, {5)g should 
converge to —1. In other words, the maximal dimensionless 
expansion of voids is —3/2 (see Bernardeau et al. 1997 for 
details). Note that the simplified formula (^ naturally sat- 
isfies this constraint. This formula is approximately valid in 
the limit ag 0. it includes contributions from all orders in 



PT. We will therefore modify it, phenomenologically incor- 
porating the (Tg-dependence. 

The formula we have found to be quite accurate for the 
whole range of values of 5 (see Figs. is 



W9 = /?K) 



i + - 

7 



where /3 is such that 



(18) 



(19) 



so it depends on the variance. This formula is relatively 
simple because the (jg-dependence factors from the 9- 
dependence. For the fit to be accurate also in the high- 
density tail, we have slightly modified the exponent [7 = 3/2 
in the expression (^)] to 



7= 

where 

e~0.3. 



^ + o /(") 



(20) 



(21) 



The normalization constraint leads to a value of /3 given 
approximately by 

1- [0.20 + 0.06/(0, A)] o-|, (22) 



/3 



for (jg below 0.8. The /(il, A) dependence exhibited here re- 
flects the slight dependence of 02 and r2 on Q that is found 
numerically. The effective coefficient a2 can be perturba- 
tively identified from the relation /3 ~ 1 — a2cr|; it varies 
from 02 = 0.23 for the low SI models to 02 = 0.26 for the 
Einstein-de Sitter Universe. 

This description reveals very accurate in describing the 
overall shape of the joint PDF in various regimes, for dif- 
ferent values of the variance and different cosmologies as 
demonstrated in Figs. |5| and ^ respectively. 

Equation ( [l7| ) together with ( p^ and (|l^) constitutes 
the main result of this paper. 



4 TOWARDS A SEPARATE 

DETERMINATION OF THE BIAS AND D. 

The description of the joint distribution of S and 9 discussed 
in the last Section has been derived for the matter density 
distribution assuming the matter density field is directly ac- 
cessible to observations. In case of biases between the galaxy 
field and the matter density field the quantitative results we 
present are not valid for the (;a/asy-density versus velocity- 
divergence relation. However, most of the qualitative fea- 
tures of the approach we propose here is expected to remain 
correct. In this Section we consider some of the new quanti- 
tative features arising if a matter-galaxy bias is present. 

Let us assume that it is possible to expand the local 
galaxy density field, 5^^'' in terms of the initial matter den- 
sity field. 



(23) 



We assume then that, from a PT point of view, the relation 
between 5^^^' and S is similar to that between S and 9. As 
a consequence, the linear bias, that is the bias parameter at 
large scale, is given by. 
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or equivalently by, 



l,gal. 



(24) 



(25) 



This identification is possible since the stochasticity is as- 
sumed to be negligible at linear order that is, 

5f'-(x) =6^^'<5i(x), (26) 

where 6^^' is a mere number. This assumption seems rather 
abrupt, but we think that it is actually completely natural. 
Indeed it does not mean at all that the final stochasticity will 
vanish. It simply means that the so called r parameter (r = 
(S'^''^- 2)1/2 ^^2^1/2) ig y^n^f-y at very large scale when 

only the linear order plays a role. This is indeed observed 
in the numerical experiments of Blanton et al. (1999). The 
second order term in Eq. ( p^ however is not assumed to be 
simply proportional to the square of the local linear density 
contrast but may contain non-local terms (as for the matter 
density contrast) or intrinsic stochasticity. Our assumption 
is even weaker than what is obtained in explicit models of 
biasing like BBKS for which the r parameter is arbitrarily 
close to unity|^ (Mo & White 1996). 

For the fields and 9 (here we reintroduce the ob- 
servable velocity divergence field 6) the forward constrained 
average reads, 



(62) 



(612)2 



2^3 



+ ■ 



(27) 



It follows from the writing of the first mom ent of the joint 
distribution function of and 9 (Eq. A13) which implies. 



(S'^'-^-e^) d^ 



2 d6l2 

with, using the Edgeworth expansion|j||. 



p{e) + . 



p{e) 



(28) 



p{e) = 



1 



211(70 

1 + 



exp 

/a'i\ 



6a 



+ . 



(29) 



Each ensemble average term appearing in this expres- 
sion can be written in terms of /(fi. A), 6^^'', Ta and observ- 
able quantities. The second moments read, 

i) = fe^^'-/(J7)(el), (30) 

(31) 



{6^) = fmiel). 

The third moment is given by, 

(e^) = -f{^i)f;{el)\ 



(32) 



^ It is unity in the limit the size of the objects is much smaller 
than the scale at which the correlation coefficients are calculated. 
II This expression is valid only when the distributions of the 
two quantities are weakly non- Gaussian, i.e., if {S^'^^-^ dt) ~ 

2{p + q-l) 



where T3 is a skewness parameter of the scaled divergence. 
This is a known quantity from Perturbation Theory that is 
independent of the cosmological parameters (with a slight 
dependence on the power spectrum index), 

96 

r3 = y-(n + 3). (33) 

In a perturbative expansion we have (iJ^^'' 9^) ~ 
((Sf el) + 2((Sf 6162) and replacing respectively Oi by 
-/(n)/6s=»i ,Sf'- and 5f'- by b^'^^ h in this expression, we 
get 



+ 



(34) 



where 5*1^' ' is the large-scale limit of the galaxy skewness. 



— 



gal. '^v 

(5g~')2' 



(35) 



This quantity can be measured directly in galaxy catalogues. 
The forward relation then reads. 



^ /(fi) 

7^ 



6 



(36) 



and the inverse relation can be easily obtained either from 
direct calculations or straightforwardly from the forward re- 
lation, 



1 ffm ^gai. fm 



6\b' 



!al. "^3 



fogal. 



T?j (5^^'-' -{5^^' ')) + ...(37) 



We stress that these relations are exact irrespective of 
bias models (linearity, scale dependence and stochasticity, 
with the restriction we mention for the linear bias). 

As can be seen in Eqs. (p6|-|3^) the coefficient in the 
quadratic terms (either a2 or 7-2) depends on b^^^' and f{fl) 
in a combination which is different from the usual one, 
/(f2)/6^^''. This is therefore a potential way of removing 
the degeneracy between the f{0,) and 6^^'' parameters. 

The aim of this paper is not to demonstrate the viability 
of such a method with real data set. This is left for future 
studies (Chodorowski, in preparation). 



5 CONCLUSION 

We have presented the results of Perturbation Theory on 
the joint relation between the density and the velocity di- 
vergence for a top-hat filter. These results extend those ob- 
tained previously for a Gaussian filter (CL97, CLPN). We 
have obtained all quantities that can be derived in third 
order perturbative calculations. That includes 

• next to leading order terms (containing so called loop 
terms) for the linear coefficients ai and ri in eqs. (^, ^); 

• the coefficients a2, 03, r2 and rs describing the con- 
strained averages expanded as third order polynomials; 
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• the coefficients &o and sq describing the constrained 
scatter up to the order four in a. These coefficients also 
contain loop terms. 

We have then used a series of numerical simulations to 
check the qualitative and quantitative behaviours predicted 
by these results. We have found a good agreement with all 
predicted qualitative behaviours. All PT results that are 
strictly at leading order (no loop terms included) have been 
found to be quantitatively well reproduced. However all the 
predictions involving one-loop calculations (that includes ai, 
ri, 60 and so coefficients) have been found to overestimate 
the results. This is particularly obvious for the linear coef- 
ficients (in Fig. ^) where the departure from the horizontal 
lines (linear predictions) is much weaker than the one loop 
order PT predictions (thin dot-dashed lines). 



France-Pologne' of CNRS/PAN. R.S. is supported by NASA 
AISRP grant NAG-3941. 



Th3 second part of our work consisted in building a 
complete description ot a joint distribution of the density 
and local velocity divergence. It appears, as we pointed out, 
that it is not possible to build the joint PDF from a proper 
Edgeworth expansion, as it should be natural in the regime 
we are. We have thus been forced to make further assump- 
tion on the statistical behaviour of some reduced quantities, 
namely that the constrained density is Gaussian distributed. 

To achieve the joint PDF construction we have taken 
advantage of the PT results that revealed accurate but have 
also extended them, using the scaling relations they sug- 
gested, with numerical fits when necessary. This is important 
in particular when one wants to deal with rare event tails 
(which are not a priori properly described by PT results). 

Our final description includes nonlinearities (such as the 
skewness) in both the one-point PDF of each quantity and 
in the expressions of the constrained means. The nonlineari- 
ties, because of the nonlocal effects they contain, also induce 
a stochasticity in this relation. We have described here its 
generic behaviour in terms of the constrained scatter and 
eventually proposed an approximate expression for the joint 
PDF providing a complete statistical description of the re- 
lation between the local divergence and the local density in 
the discussed regime. Though for technical reasons we have 
focused our presentation on top-hat filtered fields, we expect 
no significant changes for other window functions. 

Presented results may facilitate more advanced compar- 
isons of the local density-velocity data of the present day 
survey or the survey to be available in the coming years and 
provide us with a clue towards a separate determination of 
biases and Q.. 

Using as a specific example the relation between matter 
and galaxy fields, we have also discussed how local scalar 
fields employed in cosmology might be related at scales of 
a cosmological interest. In particular, we have shown that 
nonlinearities are bound to induce not only a nonlinear bias 
but a significant amount of stochasticity in their respective 
relation even if that is not present at the linear order. 
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APPENDIX A: THE JOINT 
DENSITY-DIVERGENCE PDF 
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Al General expression of the joint PDF 

The joint PDF of variables 5 and 6, P{5,9), is given by 
the inverse Fourier transform of its characteristic function, 
$. The characteristic function is related to the cumulant 
generating function, IC, by the equation 



^{it,it') = exp [IC{it,it')] . 



(Al) 



The cumulants, Hpq, from which fC is constructed, 



E ^(^^n^tr, (A2) 

{p,,j)7^{0,0) 

are given by the connected part of the joint moments 



(A3) 



The cumulants are a convenient measure of non-Gaussianity, 
since for a (bivariate) Gaussian distribution they all vanish 
for p + q > 3. 

The variables S and 9 have zero mean, hence kio = 



Koi = 0. Then from equations (Al) and (A2) we have 
1 



PiS,e) 



(2^)2 



exp — it 5 — it 



i(K20i^ + 2Kiltt' + K02i'^)] 



X exp 



.P+1>3 



dtdt' . 



(A4) 



Writing s — ^20^^, s' — ^Jj^t', introducing the standard 
variables 



X I 1/2 



and 



a I 1/2 



and the standard cumulants 



A 



PI 



^20 '^02 



(A5) 
(A6) 



yields 
1 



{2n) 
exp 



2 , , exp i(^s + J's ) — +2A11SS+S )] X 



E ^('^n^s'r 

.P+1>3 



dsds' . 



(A7) 



A2 Conditional probabilities and conditional 
moments 

The conditional probability, P{iJ,\v), reads, 

p{^^\u) = P{^,,u)/P{,y). (A8) 

The conditional moments are then the moments of P{fi\u). 
So the constrained average of /i is 

(A9) 



= / d/i/iP(^i,i.)/P(;.) 



It is possible to express this result in term of the cross- 
correlation coefficients. Using the general expression of the 
distribution in terms of the cumulant generating function 
we indeed have. 



1 



(2^)2 P(^) 
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dfj, ds ds' fi exp[—i{fj, s + ly s') + A(s, s')], (AlO) 
with 

00 

— A„„ „„ ^^^^^ 



p+q>2 



After some mathematics, we have, 

(m)|.= 2^r^y ds' ^{0,s') expH^s' + A(0,/)](A12) 
which, when re-expressed in terms of Pi^v), gives 



(">i-pJ^E(-')-^(;^)'^M^ 

9>1 



(A13) 



The higher order moments can be calculated in a similar 
way. 



APPENDIX B: CONSTRAINED AVERAGES IN 
PERTURBATION THEORY 

In perturbation theory the local density 5 and the local di- 
vergence 9 can be expanded in terms of the linear solution. 
We write these expansions. 



5 = 5i+52 + h + ■ 
and 

9 = h+92 + k + . 



(Bl) 



(B2) 



Statistical quantities of interest will then be build out of 
connected moments involving any order term of these ex- 



pansion such as {5i . 



'i)c- Taking advantage of this 



expansion it is easy to see that, in the weakly nonlinear 
regime, the standard cumulants, defined in equation (A6), 
obey for p + q> 2 the following scaling hierarchy (Fry 1984, 
Bernardeau 1992a), 



At, 



Sri 



v+q~2 



+ 0{a 



(B3) 



where a is the linear variance of 5 or, equivalently, of 9 
(we recall that Si = 9i). (The coefficients Spq in the above 
equation are called 'the hierarchical parameters'. They are 
joint skewness, joint kurtosis, and so on.) As a result the 
in (|A1? 



series in Ai , 



is equivalent to a Taylor expansion 
in a. The complete formulae can then be obtained using 
the Edgeworth expansion for the distribution P(j^). We em- 
phasize here that the Edgeworth expansion is not just 'a 
convenient form' of the PDF, but a direct consequence of 
the hierarchy Xp = Sp ct^~^ , obeyed by the standard cumu- 
lants in the weakly nonlinear regime. The derivation of the 
third-order Edgeworth expansion for a single variable is well- 
known (e.g., Longuet-Higgins 1963, Bernardeau & Kofman 
1995, Juszkiewicz et al. 1995), 



(B4) 



where Tp = Soq and Hn are Hermite polynomials. This for- 
mula is written here up to third order in a but could be 
generalized to any order. The conditional moments of jj, can 
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then be calculated straightforwardly. We express them in 
terms of quantities that are accessible to PT calculations. 
The first moment reads, 



E4Cr^ 

6 



(S3 



{S3 ~ Ti) Taa' 



(B5) 



where S3 is the skewness (properly normalized third order 
moment) of the density field, 



S3 



3{<S?<S2)c 



(B6) 



and Ts is the skewness of the scaled velocity divergence de- 
fined in an analogous way. Dividing by cr'' results in the inde- 
pendence of S3 and T3 on the normalization of the spectrum 
and, in the case of power-law spectra, smoothing scale. The 
remaining quantities are 



{{62 - e,f) 



and 



(B7) 



(B8) 



Expression (B5) is of course in agreement with the result of 



CL97. By equation (D5), the first term in this expression 



is equal to ru, resembling the case of purely Gaussian vari- 
ables, where = rv. However, nonlinearities of the fields 
in question induce corrections to the mean trend that con- 
tain also other terms linear in u. The constrained variance, 
(^Iw = (>^^)\.^ (M)t' Siven by. 



AK , (S3 - ^ o 

— + 18 + 



(AK (S3 - Ts)' 
V 2 9 



with, 



AK : 



(B9) 



(BIO) 



and in agreement with the leading-order formula in CLPN. 

The above formulas are given for standard variables, 
while for practical applications we would like to have them 
for the 'physical' ones, S and 6. The needed transformation 
involves the variances of 5 and 6 which are equal only at 
linear order. Taking this into account (se e C L97 for details) 
and rearranging the terms in equation (B5) yields expres- 
sion (^ for the mean density given the velocity divergence 
with the coefficients 



ai 



= l-f 



(S3-T3)r3 S4 
3 2 



a2 = 



03 = 



S3-T3 



S4 - (S3 - T3)r3 



The quantities E2 and S4 are given by 



E2 = 



{5202)0- {e^2}c + {Sl53},~ {6163), 



(Bll) 

(B12) 
(B13) 

(B14) 



and by Eq. (BJ). The quantity E2 is built of terms similar 



to the properly normalized weakly nonlinear corrections to 
the variance of the density and velocity divergence fields 
(Scoccimarro & Frieman 1996, Lokas et al. 1996). In the 
case of density, for example, we have to the leading order 



{5l) + 2(<5i53) 



(B15) 



On the other hand, E4 is composed of quantities similar to 
the properly normalized kurtosis which for the density field, 
to the leading order is (Bernardeau 1994a; Lokas et al. 1995) 

&{5\&l) + 4{5?53) 



S4 = 



(B16) 



The inverse relation ([7|), giving mean Q when & is kn own, is 



obtained immediately from equations (^, ( |B8[ ) and (Bll) 
(B14) via exchange of variables & and Q. We have thus 



ri = 

r2 = 

'■3 = 
where 

and 



l-f 

T3' 



S2 + 
S3 



6 ' 
El - {T3 



{T3 - S3) S3 



S3)S3 



El 
2 



'2<52)c - (<5|)c + {eie3)c - {SlS3}c 



El = 



S{0fe252)c - 3{5?52')c + {efe3)c - {5f53)c 



(BIT) 
(B18) 
(B19) 

(B20) 

(B21) 



As already stated, relations (g) and (Q) are not deter- 
ministic. The scatter around the mean given by these re- 
lations is up to terms of the order of given by expres- 
sions (^o|) and (^l|) respectively. [These expressions are next- 
to- leading-order extensions of formula ( |B9[ ) .] The coefficients 
&o and b2 in expression (llC 



60 = 
and 

b2 = 



{(52-^2)') {el{S2 



{el{S2 



2o-6 



92?)) {S3-T3f 



for (j^ig are 



e2)\ ^ (S3_ 



T3)^ 
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(B22) 



(B23) 



2cr6 9 
The coefficients 60 and &2 are invariants with respect to 
exchange of 5 with 6. This yields equation (^^. The linear 
variance of the density field is given by 

d^fc 



= {5i) = D\t) 



{2-kY 



P{k) W\kR) 



(B24) 



where D{t) is the linear growth factor [5\{x,t) = D{t)5-i{x)] 
and P{k) is the linear power spectrum of density fiuctua- 
tions. 



APPENDIX C: CALCULATION OF THE 
COEFFICIENTS FOR A TOP-HAT FILTER 

The purpose of this section is to calculate the values of the 
coefficients am, Tm and bm for the top-hat filter with Fourier 
representation 



W{kR) = 3j^{kR)-''^'j3/2{kR) , 



(CI) 
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where J3/2 is a Bessel function 



■^3/2(3;) 



(C2) 



Up to third order it is sufficient to consider the linear vari- 
ance, identical for both the density and the scaled velocity 
divergence fields, a^. 

In the following we will restrict ourselves to the case of 
scale-free spectra 



P{k) = Cfc" 



-3 < n < 1. 



(C3) 



In the case of top-hat window function and scale- free spectra 
the variance is 



a = CD^{t) ■ 



9r[(n + 3)/2]r[(l -n)/2] 



(C4) 



87r3/27?"+3r(l - n/2)r[(5 - n)/2] ' 

In the calculations of the coefficients we extensively use 
the results of Bernardeau (1994a) that led to the expressions 
for skewness and kurtosis for density and velocity fields for 
top-hat smoothing. The skewness of the density and the ve- 
locity divergence are then respectively 



34 

53 = y 
and 

26 



(C5) 



(C6) 



0.190. 



(C7) 



T3 = y - (n + 3) , 

Using equations (|C5|)-(S6) we find that contrary to the case 
of Gaussian smoothing the coefficients and r2 do not de- 
pend on the spectral index, 

4 

Finding as and rs is a more demanding task. In order to cal- 
culate E4 and E4 we need to combine the results for differ- 
ent kurtosis- type terms obtained by Bernardeau (1994a) and 
in addition find the mixed term {S1S262) /a^ ■ Using equa- 
tions (A28) and (A30) of Appendix A of Bernardeau (1994a) 



f (n + 3) + A(. + 3)^ + ^,., 
f (n + 3) + A(. + 3)^ + ^F, 

A '7 K O 

-(n-^3) + -(n + 3)V-fe^ 



we have 


2312 
441 
1352 




441 




1768 


(76 


441 




682 




189 




142 




"63" ^ 


where 





d'^fc 



^(„ + 3) + i(n + 3)^ 
|(n + 3) + i(n + 3)^ 



ii,2 
3" ' 



(C8) 



(C9) 



(27r)3 

The terms containing fc^ cancel out when we calculate E4 
5536 8 



1323 7 



(n -f 3) , 



and by combining with our equations (C5)-(C6) we finally 
get the value of the coefficient 03 which is similarly to 02 
independent of the spectral index n, 

40 

-0.0101. (Cll) 



as = 



Table CI. The coefficients ai, ri and bg in the most interest- 
ing range of spectral indices for power law spectra and top-hat 
smoothing 



n 


(ai - 1)/ct2 


(ri - 1)/<t2 


bo 


-2.0 


-0.172 


+0.0244 


0.0231 


-1.9 


-0.134 


-0.0302 


0.0260 


-1.8 


-0.0850 


-0.0991 


0.0299 


-1.7 


-0.0217 


-0.190 


0.0350 


-1.6 


-f0.0643 


-0.317 


0.0421 


-1.5 


-f0.187 


-0.507 


0.0524 


-1.4 


+0.376 


-0.825 


0.0681 


-1.3 


+0.698 


-1.45 


0.0946 


-1.2 


+1.36 


-3.06 


0.148 


-1.1 


+3.37 


-10.9 


0.309 



An analogous calculation yields 

"-^-Sl^^ (c^2) 

The values ( |C7| ) and ( pll| )-( |C12[ ) are the same as in the 
case of no smoothing, calculated by CL97 and CLPN. The 
values of the coefficients ai and ri do depend on the spec- 
tral index, their behaviour is therefore similar to the case 
of Gaussian smoothing. The calculations for a top-hat fil- 
ter can be performed in a way strictly analogous to those 
presented for a Gaussian smoothing in the Appendix C of 
CL97. The results for ai and ri as functions of the s pectral 
index are given in Table CI. From equations ( |B23| ), (C5), 
(|C6|) and (CS) it follows straightfordwardly that 



0. 



Expression (B22) can be cast to the following form 
2 



bo — a — b2 



2at 



In turn, a (Eq. B7) can be expressed as 
a = — S2 — S2 , 



(C13) 



(C14) 



(C15) 



where E2 and S2 enter definitions ( Bll ) and (B17) of ai 
and ri, respectively. Therefore, all the terms determining 
the value of the coefficient 60 have been already computed. 
The resulting value of 60 as a function of the spectral index 
is presented in the last column of Table CI. 



APPENDIX D: THE EDGEWORTH 
EXPANSION FOR JOINT DISTRIBUTION 

In th e we akly nonlinear regime, the double series in equa- 
tion (A7) is a power series in a small parameter a. A nat- 
ural idea is therefore to truncate the series at some order 
p + q = n, similarly to when deriving the Edgeworth expan- 
sion of a one -point PDF. One then expands the exponent in 
integral (A7) and retains the terms up to the order of a"~^. 



These terms are subsequently integrated using the following 
lemma: 



(CIO) _L 



(1-r^) ^''^Hp,q,(p,v,r) 



{^i^-2r^i.u + v^)/(l-r')\ 



(Dl) 



3969 



Here, Hpiqi are the bivariate Hermite polynomials. 
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{-If 



7 exp 



2(1 -r2) 



(D2) 



We will demonstrate below that unfortunately this approach 
is not strictly applicable for our case, although it could re- 
main a very useful approximation. The cumulant An is the 
correlation coefficient of the fields in question. Let us write 
r = All. We have 



1 - iaa" + 0{a*) 



(D3) 



where a is defined by equation (B7). The correlation coef- 
ficient is at the leading order equal to unity, what means 
that S and 9 are strongly correlated in the weakly nonlin- 
ear regime. For the time being, let us concentrate on the 
up-to-linear terms in the expansion of the exponent in inte- 
gral (A7). We then have 



1 



27r(l-r2)V2 



H- 



exp 



n — 3 



{n-ky.kl 



du I \duj 



2(1 



+ other terms . 



(D4) 



Including only the first term of the expansion yields the bi- 
variate Gaussian distribution, while including also the terms 
with n = 3 in the double series yields the second-order bi- 
variate Edgeworth expansion (Longuet-Higgins 1963). In our 
case r is given by equation (D3), so we can write 



/(M,.,r)=exp ^ 2(1 -r2) 

= exp [-{x' + i^^)/2\ +Oia), 

where 
_ H — 1^ 



(D5) 



(D6) 



Thus, at the leading order the dependence of / on x 
and v factors out. This suggests a useful change of vari- 
ables, il-i,!^) From the chain rule we have 

(9/a^)i„ = (ai/2. 



(a 



1/2, 



^{d/dx 

1 



-a) -\d/dx)^^ and {d/dv)^^ = {d/dv)^^ 
After some algebra, this yields 



X 



E 



1 + 

n'>'A,k<.n 

-\- other terms 



27rQi/2cr 

(-ir((M-;^r-'=^'=)c 



(n - fc)!fc! q("-'=)/2(j"-'= \dx 



d 



dp 



(D7) 



Using equations (A5) and (D5), we can cast the last expres- 
sion to the form 



and Hn's are the ordinary Hermite polynomials. The tight 
correlation between 5 and 6 causes the x and v dependence of 
each corrective term to factor. This simplification, however, 
comes at a price: we cannot formally truncate the series in 
n. The leading-order contribution to the numerator of the 
coefficient Cno is of the order of a^" , hence c„o is of the 
order of unity, for all n. The numerator of c„i is similarly of 
the order of cr^", hence c„i is of the order of a. For k > 2, 
the numerator of c^k is of the order of cr^""^, hence c^k (x 
"Yhereiore, to include all corrective terms up to, say, 
quadratic in a, we can truncate the series in equation (D8) 
at fc = 4, but we cannot truncate it at any n. Why is it so? 
The variables 5 and 6 are equal at first order, so the variable 
x oc 5—0 is fully nonlinear. Thus, v is a, weakly non-Gaussian 
variable, but a; is a fully non-Gaussian one. As a result, while 
the cumulants of u form a hierarchy in a the cumulants of 
X do not and we need an infinite number of terms in the 
Edgeworth expansion to describe the distribution of x. Since 
this distribution can be derived from the joint P{5,6), the 
joint PDF cannot be formally Edgeworth-expanded (to finite 
order) either. 



P{li,p) = 



27rQi/2c 



1 + 



E 



Cnk 



n'>Z,k<.n 

-\- other terms 



(n-fc)!fc!Q("-*=)/2 



Hn-k{x)Hk{iy) 



where 



C-nk 



{{5 - ey 



(D8) 



(D9) 
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